********************************************************************************
* Modelling of representative dataset from Air Products
* Lei Zhang
* 01/13/2014
********************************************************************************

SETS
    i   "Set of Air Products' plants"       /P1*P5/
    ii  "Set of Competitor's plants"        /PC1*PC4/
    j   "Set of Demand points(markets)"     /D1*D20/
    k   "Set of products"                   /O1/
    t   "Set of time periods"               /T1*T67/;

ALIAS(t,tt);

*-------------------------------------------------------------------------------

SCALAR M    "Big M KKT"               /10000/;
SCALAR CM   "Min amount of expansion" /90/;

SCALAR UB, LB, MB, npv1, error;

parameter y1(i,j,k,t), yc1(ii,j,k,t), cc1(ii,k,t);

$include tabledata

* TABLE DE(j,k,t)        "Demand of product each market at each time period(100ton/90day)"
* TABLE DIS(i,j)         "Distance between Air Products' plants and markets(mile)"
* TABLE DISC(ii,j)       "Distance between competitor's plants and markets(mile)"
* TABLE PRICE(k,t)       "Price of products at different time periods(MM$/100ton)"
* TABLE FALPHA(i,k,t)    "Fixed cost(MM$)"
* TABLE FBETA(i,k,t)     "Variable cost of Air Products(MM$/100ton)"
* TABLE FGAMMA(i,k,t)    "Expansion cost of Air Products(MM$/100ton)"
* TABLE FGAMMAH(i,k,t)   "First 100ton higher expansion cost part of Air Products(MM$/100ton)"
* TABLE FT(i,k,t)        "Distribution cost of Air Products(MM$/100ton/mile)"
* TABLE FCALPHA(ii,k,t)  "Fixed cost of competitor(MM$)"
* TABLE FCBETA(ii,k,t)   "Variable cost of competitor(MM$/100ton)"
* TABLE FCGAMMA(ii,k,t)  "Expansion cost of competitor(MM$/100ton)"
* TABLE FCT(ii,k,t)      "Distribution cost of competitor(MM$/100ton/mile)"

FGAMMA(i,k,t) = FGAMMA(i,k,t) / 90;
FGAMMAH(i,k,t) = FGAMMAH(i,k,t) / 90;
FCGAMMA(ii,k,t) = FCGAMMA(ii,k,t) / 90;

*-------------------------------------------------------------------------------

VARIABLES
    npv             "net present value"
    cost            "cost of Air Products and competitor";

BINARY VARIABLES
    x(i,k,t)        "if expansion happen"
    w(i,k,t)        "if first expansion happen at time t"
    xc(ii,k,t)      "if competitor's expansion happen"
    u(k,t)          "big-M for disjunction";

POSITIVE VARIABLES
    y(i,j,k,t)      "Amount of product that Air Products sells to market"
    yc(ii,j,k,t)    "Amount of product that competitor sells to market"
    c(i,k,t)        "Capacity of Air products facility"
    dc(i,k,t)       "Capacity expansion of Air products"
    z(i,k,t)        "If first expansion happen (if c > 0)"
    cc(ii,k,t)      "Capacity of competitor"
    dcc(ii,k,t)     "Capacity expansion of competitor"
    ;

    c.fx('P1','O1','T1') = 72;
    c.l('P1',k,t) = 72;
    c.fx('P2','O1','T1') = 225;
    c.l('P2',k,t) = 225;
    c.fx('P3','O1','T1') = 450;
    c.l('P3',k,t) = 450;
    c.fx('P4','O1','T1') = 0;
    c.l('P4',k,t) = 0;
    c.fx('P5','O1','T1') = 0;
    c.l('P5',k,t) = 0;

    cc.fx('PC1','O1','T1') = 270;
    cc.fx('PC2','O1','T1') = 450;
    cc.fx('PC3','O1','T1') = 90;
    cc.fx('PC4','O1','T1') = 630;
    cc.l('PC1',k,t) = 270;
    cc.l('PC2',k,t) = 450;
    cc.l('PC3',k,t) = 90;
    cc.l('PC4',k,t) = 630;

    w.fx('P1',k,'T1') = 1;
    w.fx('P2',k,'T1') = 1;
    w.fx('P3',k,'T1') = 1;

*-------------------------------------------------------------------------------

EQUATIONS
    objout          "Objective function maximize Air Products NPV"
    objin           "cost of Air Products and competitor"
    ex(i,k,t)       "amount of capacity expansion"
    epd(i,k,t)      "Investment decision in Air Products capacity expansion"
    b1(i,k,t)       "When expansion exist, first expansion exist"
    b2(i,k,t)       "relation between z and w"
    b3(i,k,t)       "z(i,k,t) <= 1"
    tmd(j,k,t)      "Demand satisfaction for all markets"
    spd(i,k,t)      "Capacity and supply of Air Products"
    spdc(ii,k,t)    "Capacity and supply of competitor"
    exc(ii,k,t)     "amount of competitor's capacity expansion"
    epdc(ii,k,t)    "Investment decision in competitor's capacity expansion"
    dj11(k,t)
    dj12(k,t)
    dj13(k,t)
    dj21(k,t)
    dj22(k,t)
    dj23(k,t)
    cut
    ;

objout..            npv =e= sum((i,j,k,t), PRICE(k,t)*y(i,j,k,t))
                        - sum((i,k,t), FALPHA(i,k,t)*z(i,k,t))
                        - sum((i,k,t), FBETA(i,k,t)*c(i,k,t))
                        - sum((i,k,t), FGAMMAH(i,k,t)*w(i,k,t) + FGAMMA(i,k,t)*dc(i,k,t))
                        - sum((i,j,k,t), FT(i,k,t)*DIS(i,j)*y(i,j,k,t));
objin..             cost =e= sum((ii,k,t), FCALPHA(ii,k,t)+FCBETA(ii,k,t)*cc(ii,k,t)+FCGAMMA(ii,k,t)*dcc(ii,k,t))
                        + sum((ii,j,k,t), FCT(ii,k,t)*DISC(ii,j)*yc(ii,j,k,t))
                        + sum((i,j,k,t), FT(i,k,t)*DIS(i,j)*y(i,j,k,t));

ex(i,k,t)..         dc(i,k,t) =e= CM * x(i,k,t);
epd(i,k,t)$(ord(t)+1 le card(t))..  c(i,k,t+1) =e= c(i,k,t) + dc(i,k,t);
b1(i,k,t)..         x(i,k,t) =l= z(i,k,t);
b2(i,k,t)..         z(i,k,t) =e= sum(tt$(ord(t) ge ord(tt)), w(i,k,tt));
b3(i,k,t)..         z(i,k,t) =l= 1;

tmd(j,k,t)..        sum(i, y(i,j,k,t)) + sum(ii, yc(ii,j,k,t)) =e= DE(j,k,t);
spd(i,k,t)..        sum(j, y(i,j,k,t)) =l= c(i,k,t);
spdc(ii,k,t)..      sum(j, yc(ii,j,k,t)) =l= cc(ii,k,t);
exc(ii,k,t)..       dcc(ii,k,t) =e= CM * xc(ii,k,t);
epdc(ii,k,t)$(ord(t)+1 le card(t))..  cc(ii,k,t+1) =e= cc(ii,k,t) + dcc(ii,k,t);

dj11(k,t)..         0.95*(sum(i, c(i,k,t)) + sum(ii, cc(ii,k,t))) - sum(j, DE(j,k,t)) =g= -M * u(k,t);
dj12(k,t)..         xc('PC1',k,t) =g= -M * u(k,t);
dj13(k,t)..         xc('PC1',k,t) =l= M * u(k,t);
dj21(k,t)..         0.95*(sum(i, c(i,k,t)) + sum(ii, cc(ii,k,t))) - sum(j, DE(j,k,t)) + 0.00000001 =l= M * (1-u(k,t));
dj22(k,t)..         xc('PC1',k,t)-1 =g= -M * (1-u(k,t));
dj23(k,t)..         xc('PC1',k,t)-1 =l= M * (1-u(k,t));

cut..               npv =g= MB;

model high_point /objout,ex,epd,b1,b2,b3,tmd,spd,spdc,exc,epdc,dj11,dj12,dj13,dj21,dj22,dj23,objin/;
model follower /objin,ex,epd,b1,b2,b3,tmd,spd,spdc,exc,epdc,dj11,dj12,dj13,dj21,dj22,dj23,objout/;
model relax /objin,objout,ex,epd,b1,b2,b3,tmd,spd,spdc,exc,epdc,dj11,dj12,dj13,dj21,dj22,dj23,cut/;

*option optcr = 0;
option reslim = 100000;

*solve high_point using mip maximizing npv;
*UB = npv.l;
*
*c.fx(i,k,t) = c.l(i,k,t);
*x.fx(i,k,t) = x.l(i,k,t);
*w.fx(i,k,t) = w.l(i,k,t);
*z.fx(i,k,t) = z.l(i,k,t);
*solve follower using mip minimizing cost;
*LB = npv.l;
*
*display UB, LB;

UB = 6035.175;
LB = 906.456;

file results /results.dat/;
put results;

while(abs(UB-LB) ge 0.5,
    
    MB = 0.5*(UB+LB);    
    
    c.up(i,k,t) = +INF;
    c.lo(i,k,t) = 0;
    x.up(i,k,t) = 1;
    x.lo(i,k,t) = 0;
    w.up(i,k,t) = 1;
    w.lo(i,k,t) = 0;
    z.up(i,k,t) = 1;
    z.lo(i,k,t) = 0;
    
    c.fx('P1','O1','T1') = 72;
    c.l('P1',k,t) = 72;
    c.fx('P2','O1','T1') = 225;
    c.l('P2',k,t) = 225;
    c.fx('P3','O1','T1') = 450;
    c.l('P3',k,t) = 450;
    c.fx('P4','O1','T1') = 0;
    c.l('P4',k,t) = 0;
    c.fx('P5','O1','T1') = 0;
    c.l('P5',k,t) = 0;

    cc.fx('PC1','O1','T1') = 270;
    cc.fx('PC2','O1','T1') = 450;
    cc.fx('PC3','O1','T1') = 90;
    cc.fx('PC4','O1','T1') = 630;
    cc.l('PC1',k,t) = 270;
    cc.l('PC2',k,t) = 450;
    cc.l('PC3',k,t) = 90;
    cc.l('PC4',k,t) = 630;

    w.fx('P1',k,'T1') = 1;
    w.fx('P2',k,'T1') = 1;
    w.fx('P3',k,'T1') = 1;
    
    solve relax using mip minimizing cost;
    
    npv1 = npv.l;
    
    if(relax.modelstat eq 1 or relax.modelstat eq 2 or relax.modelstat eq 8,
        y1(i,j,k,t) = y.l(i,j,k,t);
        yc1(ii,j,k,t) = yc.l(ii,j,k,t);
        cc1(ii,k,t) = cc.l(ii,k,t);
    
        c.fx(i,k,t) = c.l(i,k,t);
        x.fx(i,k,t) = x.l(i,k,t);
        w.fx(i,k,t) = w.l(i,k,t);
        z.fx(i,k,t) = z.l(i,k,t);
        solve follower using mip minimizing cost;
    
        error = abs(sum((i,j,k,t),y1(i,j,k,t)-y.l(i,j,k,t))+sum((ii,j,k,t),yc1(ii,j,k,t)-yc.l(ii,j,k,t))+sum((ii,k,t),cc1(ii,k,t)-cc.l(ii,k,t)));
        if(error le 0.01,
            if(LB le npv1,
                LB = npv1;
            else
                LB = MB;
            );
        else
            UB = MB;
        );
    else
        UB = MB;
    );
    
    put 'UB: ', UB;
    put '  LB: ', LB;
    put '  MB: ', MB;
    put '  stst: ', relax.modelstat/;
);

put /;

loop(i,
    loop(k,
        loop(t,
            put c.l(i,k,t)/;
        );
        put /;
    );
    put /;
);


